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ABSTRACT 

In this work we analyse the operation of a biochemical 
reactor in the repeated fed-batch mode. The reaction is assumed 
to follow Haldane kinetics i.e. it is characterised by substrate 
inhibition. The feed rate of the substrate is chosen as the 
control variable. The entire duration of the operation is divided 
into different subintervals. We optimise the system performance 
by approximating the feed flow-rate in a repeated fed-batch mode 
in two ways (i) using discrete pulses (ii) using a constant 
flow-rate over different sub- intervals. In the former the 
equations lend themselves to an analytical solution. For the 
second case we use a shooting method coupled with a non-linear 
programming technique to obtain the constant flow-rates in the 
different sub-intervals. This has been analysed for two scenarios 
(i) equal duration of sub-intervals and (ii) unequal duration of 
sub-intervals. The objective here is to maximise the biomass 
production over one cycle. 

We have also studied the applicability of our method when the 
objective is to maximise the product formed in a non repeated 
fed-batch mode. The effect of converting a non-linear constraint 
to a linear constraint on the numerical convergence has also been 
investigated. 
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Chapter 1 
INTRODUCTION 

The kinetics of biochemical fermentation reactions are 
usually characterised by substrate inhibition and product 
inhibition. The fed-batch operation of such reactors results in 
their optimal performance. Here the substrate is added to the 
reactor during the course of the reaction. The biomass is usually 
present initially in the reactor. 

The optimal operation of fermentation reactors in the 
fed-batch mode has been investigated by different methods. The 
objective here is maximising either a) the biomass production or 
b) the product formed. The classical methods applied have been 
(i) the Greens theorem, and (ii) Pontyagrin’ s maximum principle. 
Weigand (1981) considered the problem of maximising cell 
productivity in the repeated fed-batch mode of operation, for the 
case of a constant yield. He posed the problem as one where we 
have to minimise the time of operation. He obtained an analytical 
(algebraic) expression which relates the total time of operation 
with biomass concentration, for an optimal filling policy. No 
upper bound on the permissible flow-rate was assumed in the study. 

San and Stephanopolous (1984) discussed the substrate feeding 
policy in a fed-batch reactor. They determined the solution using 
the maximum principle. In a later study they (San and 
Stephanopolous (1986)) incorporated the effect of a delay in the 
growth rate and determined the optimum operating conditions. 
Cazzador (1988) studied how the initial conditions of the reactor 



determine the optimal feeding policies. These authors have 
investigated the operation of the reactor in a non-repeated 
fed-batch mode 

Bonte et al(1986) and Modak, et al (1986) determined the 
optimal feeding policies for different kinetic models. They 
developed a computational scheme to obtain the feeding policy 
which maximised a) pencillin production b)cell mass production. 

Optimal feed rate policies were determined in a non 
bio-chemical engineering context by other researchers. Levin 
(1991) discussed the problem of maximising product distribution in 
a batch reactor sustaining parallel reactions. The feed was added 
to the reactor at discrete instants of time. The equations 
governing the evolution of the system here are linear. 

Biegler (1984) used an orthogonal collocation technique to 
obtain the solution to the optimal control problem. The control 
variable was determined using sequential quadratic programming 
(SQP). This method was extended to investigate systems governed 
by differential algebraic equations, and the case where the 
control profile has discontinuities (Cuthrell and Biegler (1989), 
Vasantarajan and Biegler, (1991 )) . 

Morrison and Sargent (1984) have considered the optimisation 
problem of a multistage process, Vassliadis et al (1991) have 
discussed the solution method of multistage dynamic optimisation 
problems with path constraints. The control (usually feed-rate) 
profiles are allowed to be discontinuous at switching times (i.e 
at the junctions of different stages). These problems are solved 
using SC3P. This method requires the derivative of the objective 



function with respect to control variables. The numerical method 
is very sensitive to the evaluation of these derivatives. Rosen 
and Luus (1991) have discussed three different strategies for 
computing these derivatives. 

In the context of biochemical engineering it has been 
established that the optimal policy usually consists of a singular 
arc (Weigand (1981), Cazzador (1988)). Here the feed rate varies 
continuously with time. Such a profile is usually difficult to 
implement experimentally. In the first part of this work we 
discuss the optimisation problem of maximising biomass 
productivity in a repeated fed-batch mode of a biochemical 
reactor. 

The substrate feed-rate policy over one cycle is approximated 
as (i) discrete pulses (ii) constant feed rate over different sub 
intervals. The primary motivation behind this approach is it 
enables us to determine optimal profiles, when it is not possible 
to obtain these analytically by using classical techniques like 
maximum principle etc. These latter methods are elegant in 
analysing low dimensional systems. 

In the second part of this work we have analysed the problem 
of maximising protein production in a non repeated fed-batch 
operation. The applicability of the method proposed for this 
higher dimensional system is demonstrated. The addition of the 
substrate using discrete pulses or constant flow-rates enables us 
to numerically obtain the best "practical" approximation for the 
theoretical optimal control policies. This approximation can be 
implemented experimentally. 



Chapter 2 

PROBLEM FORMULATION (Maximising Biomass Production) 

The optimisation of the reactor in a fed batch mode consists 
of the following steps. 

(i) At the start, the reactor volume is V and contains biomass 

o 

of concentration X and substrate of concentration S 

o o 

(ii) The substrate is added over the entire cycle of reactor 
operation (t^) or till the reactor volume is filled to V^. 
The reaction is allowed to proceed for t^. 

(iii) After t^ the reactor is drained and its volume reduced to 


V 

o 

(iv) The steps (ii) - (iii) are repeated for the repeated fed 
batch mode of operation. The non-repeated operation 
terminates at the end of step (ii). 

The repeated fed-batch operation can be viewed as a cyclic 
operation of period ‘t^’ . The concentrations at the end of the 
i^^’ cycle are the initial conditions of the (i + 1)^^ cycle. 

In the repeated fed-batch mode of operation we are interested 
only in the terminal behavior of the system. This is characterised 
by the state of the system after a large number of cycles. In 
this state the initial concentrations are identical to the 
terminal concentrations of a cycle. This is in contrast to the 
non-repeated fed-batch mode, where the terminal concentration may 
not be equal to the initial concentration. The most basic 

terminal state, will be periodic with period ‘t^’ and is analogous 
to the steady-state of a CSTR. 
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In each stage, the evolution of the substrate (S), biomass 
(X) and volume (V) is determined by the evolution equations 


d (XV) 
dt 


= m(S) XV 


(la) 


d (SV) - niS) XV 
-dt = Y 


(lb) 


dV 

dt 


= F 


(Ic) 


Here /i(S) represents the dependency of the kinetics on S, Y 
the constant yield, Sj, the feed concentration, F the feed rate of 
substrate. 

A common assumption made in most of the earlier works is that 
the initial concentrations in the reactor X , S satisfy the 
stoichiometric relationship 

X = Y (S„- S ) (2a) 
o F o 

This relationship can be used to ensure us that X, S satisfy 


X = Y (S„ - S) (2b) 

r 

for all instants of time. This enables us to eliminate S in 
favour of X and reduce the order of the state , equations by one. 
This mathematical simplification is possible only when (2a) is 
satisfied for the non-repeated fed-batch operation (NRFB). Here 
the fed-batch operation is terminated at the end of step (ii). 
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For the repeated fed-batch (RFB), the system evolves to a state, 
where the relationship (2b) is true for the cycle, for 
sufficiently large N. The condition at the begining of the first 
cycle can hence be allowed to be non stoichiometric for this 
operation. We prove this in Appendix-l. For the repeated 
fed-batch mode of operation consequently it is not necessary to 
assume that the initial concentration satisfy (2a) , as was assumed 
by Weigand (1981). The concentrations X, S in the cycle 
satisfy the stoichiometric relationships (2a) for the case of 
repeated fed-batch operation for N sufficiently large. 

Eliminating the dependency on S in the kinetic expression using 
(2b) we can rewrite system (1) as 


dt 


= X m(S(X)) - 


(3a) 



(3b) 


We solve this system of equations subject to the constraints, 


0 :£ F ^ 


F 

max 


and 


V < V :£ V. 
o f 


Our objective is to determine the optimal feed-rate policy 
F(t), which maximises the productivity over a cycle, in the 
terminal state of operation. This means that we want to 



Maximise P = 


(4a) 



where is the duration of a cycle. 

Weigand (1981) obtained the optimum solution by posing this 
problem as a minimisation of time problem for a given X^. In this 
work we fix the duration of each cycle t^, and consider the 
optimisation problem as one of maximising X^. Consequently our 
problem is a fixed time problem. 

This reduces the problem to maximising P, 

P=X^ (4c) 

for a fixed V , V_, t.. 

o F I 



Chapter 3 

METHODS OF SOLUTION (MAXIMUM PRINCIPLE ) 

We now formulate the solution to this optimisation problem using 
the classical Pontyagrins maximum principle. 

The Hamiltonian H for this case is defined as 


H = 


A 

X 


Xm(S(X)) - 


F X ■ 
V 


+ A F 

V 


r \ ^ ^ 


(5) 


The evolution of the adjoint variables is governed by 


A 


X 


-A 

X 


d 

dX 


(Xp(X) - 



A 


V 



FX 

V 


subject to 


A (t.) = 1.0 

X f 


(6a) 

(6b) 


(7a) 


A (t^) = 0.0 

V f 


(7b) 


Since, the Hamiltonian is linear in the control variable F, 
the control is determined as 




^ 1 



F = F 

max 

when 

ar 





dH 


(8) 

F = F . (=0), 

min 

when 

ar 



When Hp vanishes over 

a finite 

time interval (a. 

b ) we have 


the case of singular control. 



This can occur when 


° (9) 

The control variable F can be determined in this interval 
using the information from the higher derivatives. Thus setting 

Hp = 0 (10) 

we obtain F = p.(S(X))V, For the feed rates given by this clearly, 

X = 0 

X(t) = X(a) 

The value of X, where = 0, is obtained from 

^ = 0 ( 11 ) 

A detailed derivation of the various expressions can be found 
in Appendix-2. (11) implies that the kinetic dependency of p on 
X,S must exhibit an extremum. An example of a system exhibiting 
this behaviour is when the reaction is governed by Haldane 
kinetics i.e., when we have substrate inhibition. Here 

J,(S) = (12) 

(a + I3S + yS ) 

We consider the solution to the problem for the situation 

when there is no upper bound on F i.e. (F = oo). Here it is 

permissible to obtain as large a flow-rate as desired. The 

optimal policy determined now implies we fill the reactor upto 

such that now X becomes X.from X (=X_). The reactor is then 

1 of 

filled upto using the singular arc profile along which F varies 
continuously. Along the singular arc the conditions are such that 
the reaction rate is a maximum. This reaction is then allowed to 
run in batch-mode if necessary till the end of the cycle Ct^). 



Details of this analysis can be found in Weigand (1981). 

The optimal policy of operation obtained using the maximum 
principle consists of three parts: 

(i) Here the manipulated variable F is constant at its upper 
bound 

(ii) Here the manipulated variable F is held constant at its 
lower bound. 

(iii) Here the manipulated variable varies continuously. This 
region is called the singular arc. 

It is usually impractical to vary the manipulated variable 
such as F(t) continuously as a function of time. In this work we 
discuss different possible ways in which the system can be 
operated experimentally by dividing the entire period of operation 
i.e. cycle into different interval or stages. In the i^^ stage, 
the time t is such that < t < t^ and the manipulated variable 

is assumed to be maintained constant. The reactor volume and the 
feed-rate policy in the i^^ interval are devoted V^, F^ 
respectively (Fig.l). 
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^2 ^3 
Vi V2 V3 
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Fi 

Vi 


fiN 

Vf 


^i-1 <i 


Fig.l. The variables in the different stages 

of a cycle. 
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Chapter 4 

DIFFERENT APPROXIMATIONS 

The optimal profile as determined analytically (Weigand 
(1981)) consists of three parts (i) rapid fill to an intermediate 
level (ii) singular control and (iii) batch mode operation. In 
solving our problem we determine the optimal profile by 
considering two modes of operation 


(a) We allow 

for 

the first fill 

to be a 

rapid 

fill 

and 

approximate 

the 

singular arc. 

This is 

called 

the 

ASP 


problem. 

(b) We solve for the optimisation problem by approximating the 
entire cyclic operation in discrete steps. This is called 
the total optimisation (TO) problem. 

4.1 Instantaneous Pulse Feed: 

4.1.1 The TO Problem 

We now discuss the system behaviour for the TO problem under 
different approximations. 

In this mode of operation we assume that the manipulated 

variable the feed flow rate can be infinitely large. 
Consequently, it is possible to add the substrate instantaneously 
to the reactor at various discrete time instants. This allows us 
to raise the reactor volume to through a succession of stages. 

In this mode of operation we have to determine the (i) voliame 

changes at every pulse feed and (ii) the time instants of 

addition. In Fig. 1 we depict the situation where we have N-pulses 
to raise the reactor volume from to V^, at t^, t^, .... 

Clearly, the problem reduces to determining the volumes V., , 
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^2- ••• the concentrations X(t^), X(tp 

such that X^ is maximised. A total of (2N— 1) variables have to be 
determined. The concentration X, evolves in each interval (t 

i’ 

+ from X(t^) to X(tj^_^^). The concentrations X at each 
junction point t^ experiences a discontinuity in this mode of 
operation. This is determined from the mass balance equations as 
V.X(t^^^) = V.^^XCt^^^) for i = 1,....N-1 (13a) 

Since in each interval or stage the system behaves as a batch 
reactor, the time interval of each sub-interval can be obtained as 


t. - t. 
1 + 1 1 


t . 


i+1 


dX 


X 




X)i(x) 


(13b) 


So, the total cycle time t^ is given by 


X t 


t. = 


X 


X. 


r IJ 

• 1 r 


dX ^ 




Xm(X) 


Xm(X) 

' f +1 

V > 


'fti 
^ ^1 

V J 

X 


dX 


Xm(X) 


'n-1 


The variables X(t^) can be eliminated in favour of the 
variables X(t'^) using the relations (13) to yield 


"2 f +1 

3 r +1 
^xKI 

Xm(X) 

X)I(X) 


X. 


dX 


X}i(X) 


( 14 ) 



Remembering that t^ is a constant and for a maximum the 
first derivative with respect to V^’ s and X(t]^) must equal zero, 
we obtain on differentiating with respect to X(t^), 


fi ix(t^) 




M |X(tJ 


, for i = 1, . . . N-1 


(15a) 


Differentiating with respect to yields 


M 





(15b) 


The derivatives with respect to other V^’s yield 



for i=2 N-1 (15c) 


Using these relationships we can establish 


X 





(say) 


(16) 


Assuming that there can be at most two values of X for which 
the function /i(X) attains the same values, we have 



Clearly, it follows 


V 

X — 
f 


z 


(17) 



(since otherwise the first and last steps will be of 
non-existent. ) Hence we conclude 


^f WT~., ^ 18 a) 

r 

This yields 

rz 

/ \ 
t^ = N - 1 

z 

We solve this equation for z, r,X^ such that 

mCz) = (18c) 

4.1.2 The ASP Problem 

Weigand found that the optimal control policy was filling the 

reaction to V. such that the concentration became X. (where ii has 
1 1 

an extremum) using an infinite flow-rate. The flow-rate was then 

varied such that X was maintained at X. . We now assume that the 

1 

concentration in the reactor after the first pulse is X^ i.e. 

X(t'^) = X. (19a) 

o 1 

The total number of \inknowns to be obtained is now (2N-2). 

These are X(tp, X(t 2 ) X(tjJ_^). V^, ••• and X^. The 

volume is obtained from the Junction condition at t^ as 




Using the conditions for to be a maximum, we have, 




i + 1 


V. 

■ 1 


t. 




= M 


X t 


for i = 1, N-1 




V. 

1 + 


V 


i 



X 




for i=l , . . 


From these it follows that 

M(X(t^) = fiCXCt^)) = .... = fi(X(tjJ_^)) 
Since there are at most two values of X for which fx(X) 
same values, we have 



V 

— = r (say) 
n-1 


X- 

-rz -rz . f 




dX 

4 - 

dX 

X|j(X) ^ 

Xm(X) 

X^l(X) 

-^x. 

1 

z 

z 


Where 

• fx(z) = nCrz] 


and 


^f = 


Vi 


V r 


N-1 


(20a) 

. ,N (20b) 

( 21 ) 

attains the 


(22a) 

(22b) 


(22c) 



4.2 Finite Feed-Rate (The case of equal sub- intervals) 


In the earlier approximation using pulse feeds we tacitly- 
assumed that it is possible to have an infinite rate of substrate 
addition. The reactor volume can thus change instantaneously from 
one stage to another and the volume in each stage is a constant. 
We now consider the situation where in each stage the reactor is 
operated in a semi-batch mode. Here the substrate is restricted to 
be added at a constant flow-rate. The feed rates from one 
sub-interval to another are, however, allowed to be different. 
Let be the flow-rate in the i^^ stage. The evolution in the 
i^^ stage is governed by 


HY 

= X m(S(X)) - ^ 


dt 


V. 

1 


(23a) 


dV. 


1 


dt 



(23b) 


One objective is to determine the duration of each stage ) 

and the flow rate F. in each stage such that we maximise X(t,,,) 

1 N 

which is X,. for a fixed draw down ratio (V_/V ). In this mode of 
f f o 

operation the volume of the reactor is a continuous function of 
time. Its time derivatives i.e rates of change is discontinuous 
at the junctions of the different stages. 

The concentrations X are continuous functions of time in a 
cycle. These imply the following conditions for all i 
X(t J = X(t^) 

V(t7) = V(tt) for all i (24) 

1 1 



In this mode of operation the entire operation is divided 
into N equal stages. The duration of each stage is fixed and only 
the N flow-rates have to be determined. 

4.2.1 The ASP Problem 

In this mode we increase the reactor volume from V to V 

o i 

(=V )such that the concentration X_ now becomes X. The amount of 
1 f 1 

substrate to be added can be obtained from the mass balance, 


V 


1 


V,= 

1 


^0 

X. 


(25) 


1 

This problem is the same as saying that the initial 

concentration is fixed at X. . Here however the initial volume V. 

1 1 

is unknown as the optimum X^ is unknown. 

At the end of the N intervals the reactor volume has become 
V^. The flow rates in each of the sub-intervals must satisfy 
the non-linear constraint 


\ - 


V X. 
o f 

X. 


N 

y F.At. 
L 11 


(26) 


Since V^, are constant and the time of operation 

constant the objective function reduces to (4c). 

In each stage the system evolution is governed by (23) . 


is 


a 


4.2.2 The TO problem. 

In this approximation we do not assume any instantaneous 
addition of substrate at t^ to bring the cell concentration to X^. 

So here the initial cell concentration is the same as the 
final cell concentration. The different flow-rates are to be 


determined such that we maximise X(tj^). 



Here since V^, are fixed we have the linear constraint 


N 

V - = y F.At. 

f 0 L 11 

i=l 


(27) 


The optimisation problem (ASP, TO problems) issolved using 
sequential quadratic programming. This method is based on 
linearising the non-linear constraints and constructing a convex 
quadratic function from the gradients of the objective function. 


The estimation 

of 

the 

derivatives 

plays 

a vital role in 

the 

implementation 

of 

the 

algorithm. 

Rosen 

and Luus (1991) 

have 

discussed different 

methods of estimating 

these derivatives. 

We 


have obtained these derivatives using two different approaches 

In the first method a simple finite difference scheme is 
used to obtain the derivatives. Thus 




f ^ 

F.+ AF. 

- X. 

F. 

f 

1 1 
\ 1 

f 

1 


i 


(28) 


This method involves determining X('tj^) for two different F^s 

i.e. F. + AF. , F. keeping all other F.’s (J * i) fixed. The 
1 11 j 

derivative obtained using this approach is very sensitive to the 
deviation chosen i.e, AF^. This method is hence usually not 
preferred. 

In the second approach we differentiate the model equations 

analytically and directly obtain Thus the sensitivity 

i 


equations are 



ax 

dt ( ar. 

1 


F X 

-Ir (xf'tsCX) - -i- 


X 

"v TTz ^ 


av 


dt ar. 

1 


1 


These are subject to the initial conditions, 

r ax 1 


3F. 


^Vi) = 0 


(29b) 


(30a) 


av 


(t. J = 0 


ar. ' i-i 

1 


(30b) 


Integrating (29) along with (23) yield 


9X ^ ^ 
- 5 =- at t=t. 

ar. 1 


The volume V occurring in the above equations is the volume of the 

reactor prevailing at the time instant. Thus in the ith time 

interval it is V. . 

1 

_ , , . ax(^N) . , 

To obtain — , we integrate 


ar. 

1 


dt 


ax 


ar. 

1 


ax 


X4 (X) - 1 » 


V 


"j'' av 


ar. 


V2 


(31a) 


f av 1 


dt 


5Fi, 


= 0 


(31b) 


These equations are valid for ^ ^ ^ j J ^ i 

The initial conditions for this set of equations are the same 

as those at the end conditions of the previous stage. This 

procedure is used to obtain . More detail of this technique 

i 


can be found in Levis & Kramer (1985). 


The derivatives thus obtained are used in SQP to generate the 



next iterate of A shooting method technique is used to 

ensure X^= X^, for these F^’s. This is necessary since we are 
considering the repeated fed-batch mode of operation. The 
numerical algorithm is depicted clearly as a flow sheet in Fig. 2. 


4.3 Finite feed rate (case of unequal sub-intervals) 

So far we have restricted the lengths of each stage to be 

equal. The duration of each stage was hence predetermined. We 

now relax this restriction. Now when we discretise the entire 

operation into N stages we have to determine the N flow-rates 

F^, F^ and the N stage durations 

(t^-t.) (t.-t., We now have a total of 2 N variables to 

2 1 N N-1 

be determined. In this scenario we again determine the solution 
to both problems ASP, TO, 

4.3.1 The ASP Problem 

In this case as explained before we ensure that the initial 
cell concentration is X^, by adding an instantaneous charge at 
t=tQ. The initial volume is obtained from (25). 

The entire duration of the batch and the reactor volume are 
fixed. Thisgives rise to two constraints. The non-linear 
constraint 


^0 

X. 

1 


and the linear constraint 


N 


= y F.At. 

L 11 

i=l 




N 

i=l 


(32a) 


(32b) 


The evolution of the variables in each stage is given by (23). 
4.3.2 The TO problem 


Here the initial volume of the reactor is and the 

concentration is X = X^. The two constraints now are the 

o f 

nonlinear constraint 


V 


f 



N 

y F.At. 
ik ^ ^ 


(33a) 


and the linear constraint 


N 

‘f -‘o = ,1^ ''‘i 


(33b) 


Tlie optimisation method based on SQP requires the calculation 
of the derivatives of the objective function with respect to the 
control variables. The sensitivity method to evaluate the 
derivatives with respect to the flow-rates has been discussed. We 
now elaborate on estimating the derivative with respect to stage 
durations. 

We differentiate (23) with respect to time to obtain, 


d r ax ] 
dt aAt. 
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fxp (X) 


ax 
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(34a) 



d 

dt 


aAt. 


(34b) 



These equations govern the evolution of in t. < t < t 

5At^ 1 i+1 


They are integrated subject to the initial condition 

( F. X 

= Xm (X) - 

t=tT ^ y J t=t. 


5X 


5(At. : 
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(35a) 


5V 

a(At. ) 

1 


= F. 

t=tt ^ 
1 


(35b) 


th. 

The integration of these equations upto the N stage gives 

dX (t^) 

US the required derivative > where At^ is the stage 


duration of the i^^ stage. 


4.4 Transforming a Non-Linear Constraint to a Linear Constraint: 

For a fixed final volume of reactor and a fixed time of 

operation, the stage durations and the flow-rates during the 

stages are constrained by (32a) and (33a). 

For the case where the sub-intervals i.e. stage durations 

have to be determined this is a non-linear constraint, 

This is true when the variables to be obtained are the F.*s, 

1 

At. 's. 

1 

An alternative method is to look at the total amount added 

during a stage as the process variable. Since the flow-rate in 

each stage F^ is a constant, we have 

U.= F.At. (36a) 
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The non-linear constraint now is modified to the linear 


constraint 



(36b) 


IN 


= V. « V 
f o 


The 2N variables to be obtained now are the N U.’s and the N 

1 

At^'s. The method of SQP is based on linearising non-linear 

constraints. Since the constraint in terms of the U.’s is linear 

1 

to begin with this transformation renders the numerical algorithm 
more efficient. 

The derivatives with respect to the U^’s needs to be 
obtained. This can be obtained in terms of the derivatives with 
respect to At^* s using the chain rule of differentiation. 

Let the charge added over the interval At^ be denoted by U^. 
Since the control variables are At^’s, the derivative of 

interest is evaluated for a constant value of all the other 
parameters. 

Treating as a function of At^'s, we obtain 
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(37a) 


Treating as a function of U^’s, At^’s, we obtain 
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Since is dependent upon the F^’ s and At^’s we again have 
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(38) 


Substituting this in C37b)and equating the coefficients of dAt^'s 

and dF!s we obtain 
1 



(39a) 
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(when U.’s, At.’s are control variables) 

1 1 

The derivatives required are computed using the derivatives 

evaluated with F.’s and At.’s as process variables from (29) and 
1 1 

(34). 

The simulation of the system with these linear constraints is 
such that the upper bound specified now is on U^’s, as opposed to 
F’ s in the earlier approach. It is also important while using 
the linear constraints formulation to impose a non-zero minimum 
bound on the At^’s. This is necessary as we want to avoid a time 
interval from going to zero. It is hence not possible to obtain a 
direct comparison and validation of the results using the two 


approaches. 
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Chapter 5 

MAXIMISING PROTEIN FORMATION IN A NON-REPEATED FED-BATCH 

So far we have discussed the optimum feeding strategy which 
results in maximization of biomass formed in a biochemical 
reactor. The time of operation, the initial and final reactor 
volumes were fixed. The reactor was operated in a repeated 
fed-batch mode. 

A second class of problems concerns us with the optimisation 
of the product formed in a biochemical reactor. Park and Ramirez 
(1988) considered the problem of maximizing a heterologous protein 
secreted in a fed-batch reactor, using a model host yeast SEY2102. 
The model secretory protein is SUC2-S2. 

The system evolution is governed by the set of five ordinary 
differential equations 

= A(S)(Pt-P„) - ^ P„ (40a) 

= B(S)X - ^ Pj (40b) 

X = C(S)X - |x (40c) 

S = - YC(S)X + 5(S -S) (40d) 

V F 

V =F (40e) 


with A(S) 


4.75C(S) 

.12+C(S) 


B(S) 


0.1+ s 



CCS) 


21. 87 

(S+0.4)(S+62.5) 


In these equations ( P.^) represents the level of secreted 
(total) protein in the reactor. 

The objective function is to maximize the total secreted 
protein SUC2-S2, i.e. 

J = Pj^(t^)V(t^) ( 41 ) 

Park and Ramirez (1988) solved this problem using the maximum 
principle. They considered the operation of a single cycle of the 
batch operation i.e. operation in the non-repeated fed-batch mode. 
They obtained the existence of multiple singular arcs as their 
optimal feeding strategy. They investigated the reactor 
performance for two different cycle times of 7.5 hrs and 15 hrs. 

The final volume of the reactor in this mode is not fixed and 
this results in the objective fimction being defined as in (41). 
This study was based on system performance on a single cycle of 
operation. Consequently the initial conditions of the reactor 
affect the performance. 

We extend and apply the numerical methods of the earlier 
section to analyse this problem. The algorithm discussed in 
theearlier (Fig.2)section is modified. The ‘shooting method’ step 
is deleted, since now the terminal conditions i.e. concentrations 
at the end of the cycle do not necessarily have to be equal to the 
initial conditions. We have solved this problem only for the case 
of ‘N’ equal sub-intervals. Here consequently the control 
variables to be determined are the ‘N’ F^s. Since the final 
volume is not fixed, the non-linear constraint (33a) is not valid 


29 


anymore. 

The derivative of the objective function with respect to the 

flow-rates F’ s were calculated using the sensitivity equations as 
discussed earlier. 

An important feature of this problem is at the theoretical 
analytical profile consists of two singular arcs. More 
specifically the system starts with a singular control, is 
followed by an interval where the reaction is in batch mode and is 
then again followed by a singular arc. 



Chapter 6 


RESULTS AND DISCUSSION 

We now discuss the results obtained using the different 
approximations presented in the earlier section for the two 
problems . 

The optimisation problem as posed in (3) and (4)can be solved 
using the maximum principles. Weigand (1981) obtained the 
relationship 


t 


f 


1 

M(X ) 


In 




V X- 
o f 



dX 

X m(X) 


(41) 


This relationship is derived in Appendix 2. This relates the 

time of operation t^ to X^. for a fixed V^, V at the optimum 

f f f o 

conditions. Weigand (1981) obtained this relationship for a fixed 

X^ and minimising the time t^. This relationship is also true for 

maximising X^ for a fixed t^. The variation of the control 

parameter F(t) versus time and the corresponding change in volume 

are depicted in Figures 3 a,b,. The biomass concentration X^ = 

1.95. This is the same as the initial concentration X . At time 

o 

t a pulse is added instantaneously bringing the volume upto V. 

0 

and raaintainng the concentration X at X^, where the rate is a 

maximum. The system proceeds along the singular arc for 7.922 hr 

and the final batch operation proceeds upto 10.2965 hrs.The 

parameters chosen for the simulation are a = 25, (3 = 62.9, y = 

1.0, fx = 21.87, S_ = 10, Y = 1/5, V =1.0 , V =21.0, t 

'^max F o 1 f 

=10.2965. The optimal solution X^ = 1.95. 



The optimal control policy is hence such that X, is 

maintained at where the reaction rate is a maximum. The 

singular arc is terminated once the volume of the reactor attains 

the final value V^. The reactor volume increases in the singular 

arc region till it attains V^. The reactor volume is then 

invariant till the end of the batch operation. 

In Fig. 4a we present results of the approximation of 

singular arc using discrete pulses. Here we assume that a first 

pulse is added at t such that it brings the concentration to X 

o i 

The results for N = 6 pulses, and N=16 pulses are depicted in 
Figures 4 a,b. Here we have plotted the volume of the reactor 
after the addition of each pulse. The volume changes 
discontinuously at the points where the pulses were added. These 
are the best possible approximation of Fig. 3b for agiven N. 

The duration of the stages here when the reactor volume is 
V^, ^ 3 * • ■ . equal. The first and last stages are of 

unequal duration as can be seen in Fig. 4a. The duration of the 
last batch stage compares favourably with the prediction in Fig. 
3a. 

The ratio of the volumes V. ,/V. is constant for N = 16. The 

1+1 1 

addition of the pulses is such that the biomass concentration 

after each pulse attains a constant value, X . This value X 

c c 

tends to X. as we increase N. The profile shown in Figures 4, 

1 * 

approach the analytical profile of Fig. 3b as we increase N. The 
objective function was computed as X^ = 1.9072, and 1.9142 when N 
=6, 16 respectively. In Figure 5a, b we have depicted the results 
for the total optimisation problem for the two cases of N = 6, 16. 


Here the addition of the first pulse at t^ is not constrained to 
convert the biomass from to X^. Here the duration of the first 
N-1, stages are equal. The volume ratios V^/V^ ^ are again 
constant for N = 16. The final batch operation duration matches 
the results depicted in Fig. 3b. Once again the biomass 
concentrations- after the addition of each pulse is a constant. 

This constant tends to X^ as N The biomass concentration at 

t^ for N=6, 16 ar X^ = 1.909, 1.915 respectively. These are both 
marginally better than the results of the ASP problem depicted in 
Fig. 4a, b. This is to be expected since in the TO problem we do 
not constrain the first pulse to be such that it converts X^ to 
X^. Hence the ASP problem can be viewed as the TO problem with a 
constraint. The objective function for the ASP is hence lesser 
than that of TO for a given N.In Appendix 3 we show how the 
relationship (41) can be obtained as a limiting case of (18a). 

We now discuss the results when the entire batch operation is 

divided into N stages. The stages are of equal duration and the 

flow-rates in each stage is a constant. The flow rate F is 

constrained to be less than F (10.0). For a fair comparison we 

have simulated the ASP problem, for N = 4, stages and TO problem 

for N=5 stages. Here there is an instantaneous pulse at t^, such 

that X(t^)= X.. The flow-rates (Fig. 6) for the ASP problem 
o 1 

increase as we go through successive stages. N=l,2,3. The last 
stage is operated in batch mode with F= 0. The increasing 
flow-rate captures the trend of the singular arc depicted in 
Fig. 3a. The constraint on the flow-rate F < can be 

incorported elegantly using this approach. The results for the TO 



problem are depicted in Fig. 7. Here we have used N = 5, to take 

into account the instantaneous charge at t = 0 (t ) of the ASP 

o 

problem. Here the first stage runs with an F of (1.9). There is 
than a drop in the value of F at the second stage. The flow-rate 
behaviour captures the features of the singular arc predicted 
theoretically. The flow-rate in the first stage here is not 
necessarily equal to This is only to , be expected since we 
have fixed the duration of the first stage. 

The objective function for the TO problem (1.8882) is 
lower than the for the ASP problem (1.9148). This is because 
in the ASP problem, we have no constraint on F for the first stage 
which is an instantaneous addition. 

The results of the simulations for the case of unequal time 
intervals are depicted in Fig. 8,9. The ASP problem has been 
solved for N=4. The flow-rates F^’s in each stage shows an 
increasing trend capturing the feature of the singular arc. The 
results of the TO problem are depicted in Fig. 9 for N=5. Here 
the duration of the first stage is not fixed externally. From the 
results the optimum profile is such that we have close to an 
instantaneous charge for the first-stage. Here the stage duration 
is low and the flow-rate is the maximum permissible. The 
flow-rate in the second stage drops to a low-value and increases 
in the next two stages. The last stage again is a batch mode. The 
instantaneous filling at t=0 (subject to the constraint on F), has 
been captured by our method. Once again the objective function of 
the ASP problem (1.9162) is better than that of the TO problem 


(1.91407). 



In Fig. 10a, we have shown the results of the total 
optimisation problem, when using the linear constraints. Here the 
amount of charge added in a stage is considered the control 
parameter. Here again N was chosen as 5. The optimum profile has 
been represented in terms of s again to facilitate a comparison 
with Fig. 7b. The two profiles are almost Identical as can be 
seen. The advantage with the linear constraints is that the 
constraints are identically satisfied. This results in a faster 
convergence. The number of iterations needed to converge to the 
solution here was 31 as compared to 205 with the non-liner 
constraints. The F in the first interval is 19.95020 and the 
corresponding optimum objective function is = 1.913889. 

Fig. 10b represents the X variation for the optimal flow-rate 
profile that we depict in Fig. 10a. The optimal feed rate is 
hence such that we maintain X close to X^ till the reactor is 
filled. 

The numerical approach discussed allows us to discuss the 

optimal strategy when is lower than the maximum value on the 

singular arc. We ensure that it is high enough to satisfy the 

constraints. This enables us to study how the optimal solution 

now does not emulate the singular arc. 

Figs. 11a and 11b are for the TO problem using nonlinear 

constraints. In Fig. 11a, F is 3 and for Fig. 11b F is 4. 

max max 

Here we have chosen the < max as predicted in Fig. 3a. 

This allows us to see the effect on singular profile. When F 
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is 4, F in the 4th stage is 4 while it is 4.3137 for F = 20. 

nicLx 

For F =3. the last two stages merge with each other. For F 

max inax 



It is 


= 3,4, the optimum is 1.89881 and 1.9046 respectively. 

clear from this that as we increase F , X. increases and the 

ni3.x 1 

algorithm approaches optimal profile. 

Figs. 12a, 12b and 12c are the plots of optimal control 

policy for the case where we have converted the nonlinear 

constraint to a linear one with the upper bound on u being u = 

max 

7.0, 6.3 and 5.0 respectively. The corresponding objective 

functions are 1.9052, 1.9027 and 1.8825. 

We finally simulate the system for maximising the product 

formed in a biochemical system. We investigate the problem of 
maximisation of the protein formed as studied by Park and Ramirez. 
The optimal control profile as obtained in their study consisted 
of a singular arc, followed by a batch operation and again 
followed by a singular arc (Figs. 13a and 13b} for two different 
times of operation 7.5 hours and 15 hours. We now show how the 
methodology proposed here can be used to approximate optimal 
control profile. 

We have maximised the total amount of protein formed, for a 
fixed time of operation of 7.5 hrs and 15 hrs. For t^ = 7.5 hrs, 
we have carried out simulations when the entire duration has been 
divided into equal intervals with constant F’ s in each. The 
results for N = 5, 10 are depicted in Figs. 14a, b. We clearly see 
that there is a singular arc followed by a batch and followed by a 
second singular arc. The objective function for N=5, 10 
respectively are 3.102, 3.1233. 

From t^ = 15 hrs, we show the results of our simulations in 
Figs. 15a, b, for N= 5, 10. For N=5, the flow rate shows an 



increasing trend inthe first three intervals. It then decreases 
and shows an increasing trend again in the next two stages. The 
in between batch mode is not captured here, since it lies between 
the third and fourth stages. For N=10, however the in between 
batch mode is captured in the 7th stage. Here we can clearly see 
that the two singular arcs are obtained. The objective functions 
for N=5,10 are 29.006, 29.6315 respectively. 

In Figs. 16a, b we show the variation of the objective 
function (P^V) and (P.^V) for the optimal operation for N=10. The 
variationof these quantities is similar to that observed by Park 
and Ramirez. 

The validity of the numerical scheme for this system has been 
established. The dimension of this system is larger than the 
earlier problem (3a, b) where we maximised biomass production. We 
found our algorithm to be faster for this case than in the earlier 
problem. This arises because the approach now does not the 
shooting method to ensure that the initial and final state of the 
system are equal. 

The number of divisions N has a significant effect on the 
objective function here since the singular control occurs over a 
significant portion of the cycle time. 

The computations using our algorithms were performed with a 
variety of initial guesses. We found that our system converged to 
the same optimal profile from each initial guess. 

V^-Vo 

In our work we have assumed thatF 2: — r — 

max 

This ensures us that we can fill up the reactor to in the 
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allotted time t 
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Fig6. The approximation of singular arc using 

piece-wise constant flow rate in each stage 

for equal stage duration N=4 , F =10.0 
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Fig9 . The Total optimisation problem using piece-wise 

constant flow rate in each stage , unequal 

stage durations N=5 , F =20.0 
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FiglOa. The total optimisation problem using 

linear constraints (piece-wise constant 

flow rate in each stage, unequal stage 

durations) N=5, U =20.0, At . = 0.05 
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control variable profile 




Time ( HR ) 


Fig.lOb. The Total optimisation problem using linear 

constraints (piece-wise constant flow rate in each 

stage , unequal stage durations) N=5 , U =20 0 
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= 0.05 , variation of X. 
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Figllb The total optimisation problem using 

non-linear constraints F =3.0 

max 






F ( L/HR ) 



Time ( HR ) 


Figl2b The total optimisation problem using 

linear constraints U _^„=6.3 
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Fig. 12c. The total optimisation probl^ 
linear constraints ^ x-5-0 
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Chapter 7 


CONCLUSIONS 

The use of the maximum principle to obtain solutions to 
optimal control problems is beset with many disadvantages. The 
optimal profiles are difficult to obtain theoretically in most 
cases. This is particularly the case when we have singular 
control. This is because the convergence on the solution using the 
adjoint equations is numerically difficult because of numerical 
instability. Once obtained these are difficult to implement 
experimentally. 

We have shown how we can determine numerically the best 
strategy to maximise our objective function. This strategy can be 
obtained keeping in mind the constraints, which arise in the 
experimental implementation. Thus when we have flow-rates as the 
control variable we may be able to add instantaneously or at a 
constant value. 

The advantage of the method proposed in this scheme here is 
it allows us to easily incorporate the effect of a finite value of 
an upper and lower bound on the control variable. More specifically 
it enables us to have the upper bound lower than the maximum value 
determined on the singular arc. 

We have demonstrated how we can get the best possible 
approximations to the theoretical control profile obtained from 
Pontyagrln’s maximum principle etc. 
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Appendix-1 

In this appendix we prove that the terminal state of RFB 
satisfies the stoichiometric relationship (2a). 

Dividing (la) by Y, multiplying (Ic) by (-Sp) and adding 
these to (lb) we obtain 


dt 



+ SV - VS„) 
r 


0 


This equation is valid for each cycle . 

The i*"^ cycle extends from t^ to t^ (Fig.l). To allow for 

o f 

discontinuities in the variables, the initial condition of X for 
the i^^ cycle is denoted as X(t^'*’) and the terminal condition as 
X(tp"). 

Clearly, then we have 




continuous at the different stages. The terminal state at t^ can 
be related to the initial state at using this recursive 
relationship as 



t=t 


n+ 


o 


n 


r 




+ (l-r) (l+r+r^. . +r’^ bs_ 

t=t ^ 

o 


Clearly, for r < 1, as n 


00 , 


n 


r 


^ 0 and we obtain 


+ S 


t=t 


n+ 


The initial concentrations in the terminal state of the RFB hence 


satisfy the stoichiometric relationship 2a. 
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Appendix-2 

In this appendix we present the derivation of the singular 

arc. 

Hp — 0 for a finite time interval,, say some a t i b and 
we cannot determine the value of F from = 0. we will have 
singular control. To determine this, we get higher order time 
derivative of and solve them simultaneously. In this case 

Ay M' (X) 

TT _ A 


iij, = A^ X^[-2F Ii' (X) + p' (X) V p' (X) 

- VX fi' (X) - FX (H"(X) 

+ p(X) VX |i"(x)] / 

On reduction of above two expressions, we get 
Aj^ 0, m"(X) ^ 0, M' (X) = 0 
F = m(X)V and A^ - (X/V) A^ = 0 

So we will have singular control only when p' (X) =0. If we 
substitute this value of F in the state equation we get 
X = 0 
dV 

■TT = m(X. )V (where X. is the solution of u' (X) = 0). 
dt 1 1 

V = V(a) exp (ji(X^) ) (t-a) 

F = |i(X^) V(a) exp. (pCX^) ) (t-a) 

Thus our total time of operation is divided into two parts: 

(i) Time in filling the reactor along the singular arc 

(ii) Time for batch operation 
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and is given by 

^fill ^batch 
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Appendix-3 

In this appendix we show how the optimal control policy for 

N-pulses converges to the theoretical optimal policy for N xo. 

Expression for in total optimisation problem for (N+1) pulses 
is given by 
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Using L-Ho’pital rule 
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This is identical to the expression of Weigand as derived in 
Appendix-2. 

Now for t^ to be minimum for a fixed X^ 



= 0 


This yields 


where, 
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in terms of the kinetic parameters in C12). 



